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^ ^ , Abstract 

' In this paper we consider systems of partial (multidimensional) linear difference equations. Specif- 

, ically, such systems arise in scientific computing under discretization of linear partial differential 

equations and in computational high energy physics as recurrence relations for multiloop Feynman 
integrals. The most universal algorithmic tool for investigation of linear difference systems is based 
\^ , on their transformation into an equivalent Grobner basis form. We present an algorithm for this 

^0 ' transformation implemented in Maple. The algorithm and its implementation can be applied to au- 

I tomatic generation of difference schemes for linear partial differential equations and to reduction of 

O ■ Feynman integrals. Some illustrative examples are given. 

1 Introduction 

t:;;^- I Let N>o and N>o be the sets of positive and nonnegative integers. Q be the set of rational numbers, 

' Y := { (xi, . . . , Xn)} | j = 1, . . . , m m,n N>o} be the set of functions in n-variables, and 9i be the 

right-shift operator for the i-th variable: 9i o y{xi, . . . , Xn) vi^i, . . . , Xi + 1, . . . , Xn) ■ For the power 
_ products 9f^^ ■ ■ ■ 9!^" of the shift operators we shall use the multiindex notation 0^ where yu := {/ii, . . . , /i„} 

«) I {fi e N"q) with fi -=^2 Ms- The set of all such operator products will be denoted by 6. 

""-^ ■ And then the most general form of a system oi K ^ N>o partial (n > 1) and multivariate (m > 1) 

^ , linear difference equations is given by 



ao + J2J2''^r,-^k°y' k = l,...,K, K&Q, (1) 

H ; 

' where all sums are finite and coefficients oo, akj-. v may depend on the variables X := {.ti, . . . , Xn\ and on 

a finite set of parameters C := {ci, . . .}. Hereafter we shall assume that all coefficients in ([1]) are rational 
functions of the variables and parameters with integer coefficients: 

ao,afc,;. e Q(^UC). (2) 

This restriction on the coefficients allows to apply algorithmic technique of the next section. 

It is well-known that, except very simple cases, systems of form (IT][2]) do not admit exact solutions and 
rather weakly studied in the literature [1]. However, such systems play a fundamental role in a number 
of important applications for instance in: 

Scientific Computing: Numerical solving of linear partial differential systems [2] with rational function 
coefficients. Recently, it was observed [3] that one can automatically generate finite-difference 
schemes for such partial differential equations (PDEs) by eliminating partial derivatives from certain 
linear partial and multivariate difference systems. In so doing, for homogeneous PDEs whose 
coefficients may also be rational functions of parameters one deals with systems of form with 
flo = 0. 
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Computational High Energy Physics: Reduction of multiloop Feynman integrals [4]. These inte- 
grals, after a proper fixed right shift of the variables in X satisfy the univariate system of partial 
difference equations (recurrence relations) [5] whose rational function coefficients depend on such 
physical parameters as the space-time dimension, masses and external momenta. The problem is 
to reduce the integrals to be evaluated to a minimal set of basic or master integrals, i.e., those 
integrals which are independent modulo the difference system, and then to express other integrals 
in terms of the basic ones. 

Computational Economics: Characterization of economic behavior in macroeconomics [6]. Here for 
some macroeconomic problems one has to solve system ([T]) with constant parametric coefficients. 

To investigate or to solve difference systems ([T]) with rational coefficients ([2]) one can use the universal 
algorithmic Grobner bases method invented about 40 years ago by Buchberger [7] for systems of multi- 
variate commutative polynomials generating polynomial equations [8] . The main idea of this method is to 
rewrite the initial system of equations into a certain equivalent form called a Grobner basis which makes 
easier investigation of the system and its solving. The underlying Buchberger's algorithm [8] built-in all 
modern general-purpose computer algebra systems such as Maple, Mathematica and others. 

On the basis of research made to date, the Grobner bases theory was extended to some "weakly" 
noncommutative polynomials as well as to linear differential or difference polynomials and operators [9l 
110] . Generally, however, the noncommutative and nonlinear differential or difference Grobner bases may 
not exist (be infinite) . For difference systems ([T][2]) Grobner base are always finite and can be constructed 
by Buchberger's algorithm straightforwardly translated to difference algebra |10j . 

Recently [3llllj we presented the difference form of our polynomial algorithm devised in [12j . improved 
in [13] and specialized to so-called Janet and Janet-like monomial divisions [14] which go back to the 
constructive ideas of French mathematician Janet [15]. The algorithm constructs a Janet (-like) basis [14] 
which is also a Grobner basis. Though generally Janet bases [T^j and Janet-like bases are redundant as 
Grobner ones, the algorithm in its improved version [S] [13] allows also to output reduced Grobner bases 
without any additional computational costs. The implementation [16j of the algorithm in Maple allows 
a user to compute linear difference Janet(-like Grobner) bases. 

In the present paper we describe briefly a simple version of the Janet division algorithm (Sect. 2) and 
consider its application to the above listed problems from scientific computing (Sect. 3) and computational 
high energy physics (Sect. 4). Both problems are purely algebraic and can be completely solved with the 
use of Grobner bases. We illustrate this fact by simple examples. Our presentation is addressed to 
non-algebraists. By this reason we slightly abuse algebraic terminology and refer to the references in 
bibliography for more precious definitions and notions. 

2 Transformation to Grobner Basis 

In this section we define the concept of a Grobner basis form for the difference system ([T][2]) and present 
an algorithm for its computation. The Grobner basis form of system ([1]) is defined by a ranking >- (linear 
order on) of terms 9^ oy^ and such that for all i, j, fc, /i, v the following holds: 

e,e^' o y^ y r o y^ , r o y^ yO^o y^ 9^0^' o y^ ^ 9,6" o y^ . 

If l/^l >- =^ 9^^ o yi >- 9'^ oy'' the ranking is called orderly, li j > k =^ 9^^ oy^ )^ 9'^ o y^ the ranking 
is called elimination. 

Denote by fk the left-hand side of the fc-th equation in ([1]) and by F := {/i, . . . the set of 

all the left-hand sides in the system. Fixing a ranking >- provides every f £ F with the leading term 
lt(/) := a-d oy^ G Q, a 0) and leading coefficient lc(/) := a. Furthermore, denote R D F the set of 
all right-hand sides ,f ^ for linear difference equations / = which are consequences of system ([T][2]). 
F is called generating set or basis of R (denotation: R =< F >). In that follows we shall assume that, 
given a ranking all / S i? are normalized, that is, divided by their leading coefficients. If C R, then 
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lt(F) will denote the set of the leading terms and It j{F) will denote its subset for function . Therefore, 

lt(i^) = U™ilt,(F). 

Now we are ready to define a Grobner basis for given F and ranking >- as a finite subset G C R =< F > 
such that R ^< G > and 

Vf eR, 3g(EG,9eQ : lt(/) ^0o \m{g) . (3) 

It follows that the leading term of every / G i? is reducible modulo G and yields the head reduction: 

f ^ r f ~0og, feR. 
g 

If /' 7^ 0, then its leading term is again reducible modulo G. And then by repeating the reduction finitely 
many times [SI [9l \W\ we obtain / — > 0. Generally, if a linear difference expression h (not necessarily 

G 

from R) contains a term u with coefficient c ^ such that u — cdo lt(/) for some e 8 and / G F C M, 
then h can be reduced: 

h ^ h' ■.= h- c-do f . (4) 
g 

By applying the reduction finitely many times, one obtains a polynomial h which is either zero or such 
that all its (nonzero) terms are irreducible modulo set F . In both cases h is said to be in the normal form 
modulo F (denotation: h ~ NF{h, F)). A Grobner basis G is called reduced if g ~ NF{g, G \ {g}) for 
every g G G. 

In our algorithmic construction of reduced Grobner bases we shall use a restricted set of reductions 
called Janet reductions (cf. [13| ) and defined as follows. 

For a finite set F and a ranking we partition every set Itfe(-F) into groups labeled hy do, ... ,di G N>o, 
(0 < i < n). Here [0]^ := Itfe(F) and for i > the group [do, . . . , di]k is defined as 

[do, . . . , d.i]k := {u G Itfc(F) \do = 0, dj = degj{u), 1 < j < i} 

where dcgj(0^ o y'^) := /ij. Operator 6i is called J (anet)-multiplicative for / G F if lt(/) G [do, . . . , di-i] 
and degj(u) = max{degj(u) | v G [dp, . . . , c?i_i]}. Denote by Mj{f,F) the set of J (anet)-multiplicative 
shift operators for / G -F, the complement set {0i, . . . ,9n} \ Mj{f,F) of J (anet)-nonmultiplicative shift 
operators by NMj{f,F) and the set of all possible power products of J— multiplicative operators (in- 
cluding identity operator) by J{f,F). It is clear that J{f,F) C 0. 
A finite set G G i? =< F > is called a Janet basis (cf.|T3]) if 

Vf eR,3geG,0eJ{g,G) : lt{f) = eo\t{g). (5) 

Similarly to a J— reduction is defined as 

h h' := h- cdo f , i}eJ{f,F), (6) 
g 

for a polynomial h G R containing a term u with coefficient c ^ satisfying u — ct?olt(/) for some f G F 
and I? G J(/,F). 

Since J— reducibility ^ implies the Grobner rcducibility ([4]), a Janet basis satisfying ([5]) is also a 
Grobner basis. The converse is generally not true, that is, not every Grobner basis is Janet one. The 
algorithmic characterization of a Janet basis G is the following condition (cf. T3]): 

Vg G G, ee NM{g, G) : NFj{0 og,G) = 0. (7) 

which is a cornerstone of the below algorithm for construction of Janet bases ([3]). 

This algorithm is a translation (with some minor modifications) of the polynomial algorithm in [17] 
into the difference case. Due to the normalization of h done at Step. 15 before insertion of h into the 
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intermediate basis G, the algorithm outputs the minimal and normalized Janet basis which is uniquely 
defined by an input difference system F and a ranking [12] . Correctness and termination of the difference 
algorithm immediately follow from those for its polynomial counterpart [121 US]- Algorithm Janet- 
Basis implemented in its improved form [3] as the Maple package LDA (abbreviates Linear Difference 
Algebra) [16j . and in the next two sections computation with the package is illustrated by examples. 



Algorithm: JanetBasis(F, )^) 



Input: F, a finite linear difference set; y, a ranking 


Output: G, a Janet basis of < > 


1 


choose f <E F with the lowest lt(/) w.r.t. >- 


2 


G {/} 


3 


Q:^F\G 


4 


do 


5 




6 


while Qt^% and /i = do 


7 


choose p <E Q with the lowest \t{p) w.r.t. >- 


8 


Q-^Q\ {p} 


9 


h:^NFj{p,G) 


10 


od 


11 


if /i 7^ then 


12 


for all g e G such that \i{g) = 6*^ o \i{h), |^| > do 


13 


Q:=gu{.9}; G := G \ {g} 


14 


od 


15 


G ■.^G\J{h/\c{h)} 


16 


Q:^Q\j{6Pog\geG, 9^ & NMj{g,G)} 


17 


fi 


18 


od while 


19 


return G 



3 Generation of Difference Schemes 

In paper [3] an algorithmic approach was developed to construct finite-difference schemes for linear PDEs 
in two independent variables and on uniform orthogonal grids with possibly distinct mesh steps for x 
and y. We outline here the main idea of the approach and refer to [3] for more details. In so doing, we 
restrict our consideration by scalar equations of order > 2 which admit the conservation law form 

dV ^dW _^ 
dx dy 

where V and W are functions of independent variables x,y^ dependent variable u(x,y) and its partial 
derivatives Ux,Uy, u^x, ■ ■ ■■ Differential equation ([8]) can be rewritten in the integral form 

-Wdx + Vdy = (9) 



r 



which is valid for arbitrary closed contour F. Discretization of ^ instead of ^ is natural for preserving 
the conservation low at the discrete level (conservative scheme). 
Denote the grid values of function u(x, y) and its derivatives by 

Ujk ■= u{^j,yk), {ux)j k ■= Uxixj,yk), {uy)j k --^ Uy{xj,yk), (uxx)ik ■= Uxx{xj,yk), ■■■ , (10) 
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Figure 1: Integration contour on grid 



and fix some integration contour F in ^ on the grid. To be specific, let us clioose the following simple 
rectangular contour 

Now we add to the integral equation ([9]) for the rectangular contour of Fig. [T] all the related (exact) 
integral relations between u{x,y) and its partial derivatives: 

f yk+2 

J Uxdx = u{xj+2-,v) - u{x-i,y) , / Uydy = u{x,yk+2) ~ u{x,yu) , 

Xj yk 

^3+2 yk+2 

J Uxxdx = Uxixj+2,y) - Ux{xj,y) , J Uxydy = Ux{x,yk+2) - Ux{x,yk) , 
Xj yk 

Our purpose is to obtain a difference scheme for Uj k from a proper discretization of integral equations 
^ and relations (fTT|) . To do that one should use as many relations in pT|) as the number of all proper 
derivatives of u up to the maximal orders of their occurrence in the integrand of Then the difference 
scheme can be obtained by an algebraic difference elimination of all discrete proper partial derivatives 
in list (fTO|) from the combined system (|9llip . The algebraic elimination can be achieved by computing 
a Grobner or Janet basis for the last system and a suitable elimination ranking (see Sect. 3) satisfying 

Ujk -< (ux)jk -< {uxx)jk -<•■■• 

Therefore, to construct an initial system of discrete equations for the following difference elimination, 
it suffices to approximate numerically the contour integral © for the chosen contour of Fig. [T] together 
with the integral relations pT|) in terms of the grid unknowns (|10p . For this purpose one can choose 
various quadrature formulas for these integrals, and the difference scheme obtained may depend on the 
choice. For simplicity sake we apply here for all the integrals in ([9]) and [TT|) the simplest rectangle 
(midpoint) rule: 

{Wj+i k+2 - Wj+ik) ■ hi + {Vj+2 fc+i - Vj fc+i) • /i.2 = , 
{ux)j+-i k ■ 2hi — Uj+2k ~ Uj k , 

{Uy)jk+-1 ■ 2h2 = Ujk+-2 ~ Ujk , 



where hi :— Xj+i — Xj and /12 yk+i — yk arc the grid mesh steps for our uniform orthogonal grid. 

For linear difference system (|12p Janet ( Grobner ) basis exists for any ranking, and, hence, the 
elimination can be performed by applying the above algorithm JanetBasis. To illustrate this algorithmic 
procedure for the difference schemes generation consider a simple example of the Heat equation in its 
conservation law form [3]: 

Ut + auxx = => Q) —auxdt + udx — . (13) 
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Figure 2: Integration contour for the Heat Equation 



where a is a symbohc parameter. The integrand in (|13p contains the only partial derivative m^. Hence, 
we need to add the only integral relation 



Now consider u(x,t) and Ux{x,t) on the uniform orthogonal grid with the spatial mesh step h and 
the temporal mesh step r, and choose the simplest contour shown in Fig. [2l As this takes place, we can 
approximate the integral of Ux{x, t) over x in p3m4|) on the grid points by the rectangular or trapezoidal 
rules. 

Then, applying the midpoint rule for the contour integral and the trapezoidal rule for the relation 
integral we find two difference equations for two dependent variables w, u^'. 



Furthermore, we show how to generate a finite-difference scheme for the Heat equation (|13p by using 
the Maple package LDA [TB] : 

> with (LDA) : 

First, we enter the independent and the dependent variables for the problem. 

> ivar := [j,k]; dvar := [ux,u] : 

Second, we translate (jisp into the input format of the main command JanetBasis in the package. 

> L:=[a*t/2*(ux(j ,k)+ux(j+l ,k)-ux(j ,k+2) -ux(j+l ,k+2) ) -2*h* (u( j+1 ,k+l) -u( j ,k+l)) , 

> h/2*(ux(j,k+l)+u(j,k)-u(j,k+l)+u(j,k))] : 

Third, we compute the (minimal) Janet basis for L w.r.t. an elimination ranking with Ux )~ uto eliminate 
the partial derivative from the system (jlSp . This ranking is chosen by using option 2 as below; in so 
doing we output only the element in Janet basis which does not contain Ux- 

> JanetBasis (L, ivar, dvar, 2) [1] [1] ; 

-2atu{j ,k + 1) + hat u(j ,k) + 2atu{j,k) + 2atu{j,k + 3)-hatu{j,k + 2)-2at u(j ,k + 2) 
+ 2atu{j + l,k + 3) - hatu{j + I, k + 2) - 2 atn{j + 1, k + 2) + Ah'^ u{j + 1, k + 2) 
-Ah'^u{j ,k + 2) -2atu{j + l,k+l) + hatn{j + l,k) + 2atn{j + l,k) 




(14) 




(15) 
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Thereby, we obtain the classical Crank-Nicolson scheme 

if in the above Maple output one shifts the second index by -1 and uses the first index as a superscript. 
The same scheme is also obtained for the midpoint integration method applied to ((TJ]). 

4 Reduction of Feynman Integrals 

Evaluation of Feynman integrals is the cornerstone step of perturbative computations in elementary 
particle physics [4]. Consider, for example, a typical scalar L— loop integral with E external legs: 

X{yi / d^ki ■ ■ ■ cf'kL ^ (16) 



which corresponds to n internal lines in the related Feynman diagram. Integration for every loop mo- 
mentum ki is done over the space-time of dimension c? = 4 — e where e is the parameter of dimensional 
regularization [TS]. The denominator Dj for the j-th internal line with mass rrij is given by Dj := Pj—m^. 
Here the line momenta pj are linearly expressed in terms of the loop momenta ki (i ~ I, . . . , L) and ex- 
ternal momenta Qs (s = I, . . . , E) a,s 

L E 

s=l t=l 

Consider the combined set of L + E vectors 

ka, a = 1, . . . , L , 
qa-L, a^L + l,...,L + E. 

Recurrence relations for integral p6p are derived by the integration-by-parts method [S] whose main idea 
is to use the integral identities (cf. [41IT9]) 

together with the vector identities 

2p.gj = (k + qjf - (p- - mj) - {q] + m^) . (18) 

Integral identities p7)) follow from an observation that any integral of d/dki{. . .) vanishes since there are 
no surface terms in dimensional regularization (cf. [20j). 

As a simple example consider one- loop propagator diagram of Fig. [3] with external momentum q and with 
one massive and another massless internal lines. This example was studied, for example, in [H 12 Ij and 
used already in [16| . The corresponding Feynman integral (|16p is given by 



^K'^2)- / T^I^TTv^, D,:={k + q)\ D^ '.^ k\ (19) 



d'^k 

In accordance to (jl7p . there are two independent integral identities. Their integrands are 



k 



q 




k + q 



q 



Figure 3: One- loop propagator diagram 



Taking into account the identity 2k ■ q ~ {k + q)^ — (fc^ — m?) — [q^ + m?) of type (jlSp and equality 
d/dk ■ k = d wc obtain the difference system 



where 6\ and 62 denote the left-shift operators for indices 1^1 and 1^2, respectively. 

Now wc construct the minimal set of master or basic integrals for the two-indexed family (fT9|) of 
Feynman integrals by applying the Maple package LDA |16| with the input denotations k := vi^n V2 
and f{k + l,n-|- 1) := l{vi,V2): 

> ivar: = [k,n]: dvar: = [f]: 

Then, we enter the recurrence relations (|20p . 

> L : = [ (d-k-2*n) *f (k+1 ,n+l) -k*f (k+2 ,n)+k* (q'2-m"2) *f (k+2 ,n+l) -2*m~2*n*f (k+1 ,n+2) , 

> (k-n)*f (k+l,n+l)+k*(q~2-m'2)*f (k+2,n+l)-k*f (k+2,n)+n*f (k,n+2)- 

> n*(q~2+m~2)*f (k+l,n+2)] : 

As the next step we compute a Janet basis for an orderly ranking (Sect. 2) induced by 9i )^ 02- 

> JB:=JanetBasis(L,ivar,dvar) : 

In order to compute the set of master integrals we have to take into account that Z(!^i,//2) = for 
1^2 < [H [20] . This extra boundary information is input as 

> AddRelation(f (k+j ,n)=0 , ivar ,dvar) : 

Master integrals are those f{k,n) which are independent modulo all the consequences R (see Sect. 2) 
of (j20|) . Thereby, the master integrals are easily determined via the leading terms of the Janet basis. 
Namely, one has to determine those /(fc,n) that are not expressible as the action of a power product 
^i^^2^ (^i,/Z2 £ N>o) on a leading term in the Janet basis (cf. definition dS])). 

The set of master integrals is computed by invoking the command: 

> ResidueClassBasis(ivar,dvar) ; 



Now any integral f{k + i,n + j) can be explicitly expressed as a linear combination of the master 
integrals whose coefficients are rational functions in parameters d,q^ and m^. The explicit expression is 
obtained algorithmically by applying the Grobner or Janet reductions described in Sect. 2. In LDA the 
Janet reductions are performed. To show the output of such an expression for f{k + 3, n -f- 2) and make 
the output more compact we let m = and show the underlying piece of the Maple code: 



( 



[rf - z^i - 2,^2 - i^iOiOi + i^i{q^ - 7712)6*1 - 2171^1^202] oI{v^, 1^2) = , 

[vi -U2 + Mq^ - m^Wi - lyioA^'^ + V20['^^e2 - 1^2(9' + ^2)^2] o2:(zyi, 1.2) = , 



(20) 



[f(fc, n + 1), f(fc, n + 2), f(fc + 1, n + 1)] 
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> m:=0: J:=JanetBasis(L,ivar,dvar) : 

> AddRelationCf (k,n+j )=0 , ivar ,dvar) : 

> ResidueClassBasis(ivar,dvar) ; 

[f(/c + l, n+1)] 

> InvReduceCf (k+3,n+2) ,J,"F") ; 

-{{d - 2 - 2 k - 27j) {d - 4: - 2 k - 2n){d - 2 - k - n) {d - 3 - k - n) {-n - 2 k + d) f{k + 1, n + 1))/ 
{q^ {k + l){-2k + d-4)k{-2k-2 + d) n) 

In the massless case (to = 0) a new extra relation f(k, n + j) equivalent to the boundary condition 
'^{vi, V2) = for vi <0 \s added that yields the only master integral f{k + 1, n + 1). The last command 
invokes the LDA procedure that computes the J— normal form of f{k + 3, n + 2) modulo the Janet basis. 
This normal form just represents /(fc + 3, n + 2) in terms of the master integral. Option "F" provides 
factorization of the numerator and denominator in the output rational function coefficient. It should be 
noted that, since integral 'I{vi^V2) is non-vanishing only when both its indices are positive, the master 
integral can be identified with X(l, 1). 



5 Conclusion 

We shown above that the Grobner bases technique can be applied to generate difference schemes for linear 
PDEs and to reduce multiloop Feynman integrals. Each of our simple illustrative examples of Sect. 3 and 
4 needs less than 1 second of computing time on an 1.7 Mhz personal computer with 512 Mb RAM. 
Larger examples, however, can require much more computer resources since complexity of a Grobner 
basis computation is at least singly exponential, and may be even doubly exponential, in a number of 
variables [9l[22]. Besides, blowing-up of intermediate coefficients, especially in the presence of parameters, 
as in the case of recurrence relations for Feynman integrals, is a serious obstacle in the practice. That is 
why to apply Grobner bases to multivariate and multiparametric problems one has not only to optimize 
and improve the underlying algorithms and data structures but also to implement them in lower level 
languages than Maple or Mathematica. Our Janet division algorithms [13] have already implemented 
in C and C-|— I- (see the Web page [13]) for commutative polynomials, and extension of these codes to 
differential and difference equations is planned for the coming years. As it is argued in [211 [24], differential 
Grobner bases can also be applied to reduction of Feynman integrals. A practical specialization of the 
Grobner bases ideas to reduction of Feynman integrals was suggested recently in [551 US] where the whole 
index space for integrals is partitioned into so-called sectors in accordance to the extra boundary 
conditions. Then in every sector a certain kind of a Grobncr-like basis is constructed. As to an extended 
discussion of generating difference schemes by means of Grobner bases we refer to our recent paper [3] . 
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